arXiv:1505.04598vl [astro-ph.GA] 18 May 2015 


Mon. Not. R. Astron. Soc. 000, ??-?? (2015) Printed Tuesday 19^^ May, 2015 (MN BTgX style file v2.2) 


Spiral density waves in the outer galactic gaseous discs 


S.A. Khoperskov^’^’^* and G. Bertin^ 

^ Dipartimento di Fisica, Universitd degli Studi di Milano, via Celoria 16, 1-20133 Milano, Italy 
“^Institute of Astronomy, Russian Academy of Sciences, Pyatnitskaya st., J^.8, 119017 Moscow, Russia 

^Sternberg Astronomical Institute, Moscow M.V. Lomonosov State University, Universitetskij pr., 13, 119992 Moscow, Russia 


Tuesday 19^^ May, 2015 


ABSTRACT 

Deep HI observations of the outer parts of disc galaxies demonstrate the frequent 
presence of extended, well-developed spiral arms far beyond the optical radius. To 
understand the nature and the origin of such outer spiral structure, we investigate the 
propagation in the outer gaseous disc of large-scale spiral waves excited in the bright 
optical disc. Using hydro dynamical simulations, we show that non-axisymmetric den¬ 
sity waves, penetrating in the gas through the outer Lindblad resonance, can exhibit 
relatively regular spiral structures outside the bright optical stellar disc. For low- 
amplitude structures, the results of numerical simulations match the predictions of a 
simple WKB linear theory. The amplitude of spiral structure increases rapidly with 
radius. Beyond ~ 2 optical radii, spirals become nonlinear (the linear theory becomes 
quantitatively and qualitatively inadequate) and unstable to Kelvin-Helmholtz insta¬ 
bility. In numerical simulations, in models for which gas is available very far out, spiral 
arms can extend out to 25 disc scale-lengths. A comparison between the properties 
of the models we have investigated and the observed properties of individual galaxies 
may shed light into the problem of the amount and distribution of dark matter in the 
outer halo. 
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1 INTRODUCTION 

For several galaxies, deep HI images demonstrate the pres¬ 
ence of large-scale spiral arms in extended gaseous discs 
well outside the bright stellar component (e.g., see Sancisi 


|et al. 2008). The most spectacular cases are NGC 1512/1510, 


NGC 5055, NGG 6744, NGG 6946, and NGG 5236. Some spi¬ 
rals exhibit a rather complex morphology (e.g., NGG 6946; 
Boomsma et al. 2008| ); yet, symmetric grand-design struc¬ 
tures are not so rare (e.g., NGG 1512; [Koribalski Lopez- 
Sanchez||2009 ). 

UV images of the extended galactic (XUV) gaseous disc 
point to ongoing star formation even in such outermost re¬ 
gions (Bigiel et al. 2008). In particular, GALEX data provide 


evidence for weak but signihcant star formation in such very- 
low gas-density environments. HII regions have been noted 


far away, outside the bright optical disc (Karachentsev et al. 


2011). In addition, the gaseous layer in the disc plane is char¬ 
acterized by a rather high velocity dispersion 1 — 10 km s~^. 
An interpretation of these hndings is required. 

A comparison of HI and UV brightness of the ex¬ 
tended outer disk shows that spiral arms are a necessary, 
but not sufhcient, requirement for star formation (Barnes 
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et al. 2012). Interestingly, although with only little atten¬ 


tion paid to the dynamical origin of spiral structure, the 
simulations by Bush et al. ( 2008| (see also Bush et al.|2010 ) 
show that organized compression regions and hlamentary 
structures should occur frequently, at least if the extended 
gas is characterized by a constant surface density at the level 
5—10 Mq pc“^ (which is on the high side, with respect to HI 
observations; note also that this implies a much shorter de¬ 
pletion time scale than expected ( Bigiel et al.||2010| )). How¬ 
ever, in spite of the fact that there are many examples of 
extended UV discs, we still lack hrm theoretical models of 
star formation outside the optical radius in external galaxies. 
This difhculty is largely related to our poor understanding 
of the physical conditions in those regions, both in terms of 
general properties of the interstellar medium and of the re¬ 
sulting stellar populations ( Koda et al.||2012 Barnes et al. 


2013 ). It appears that, much like in some dwarf galaxies, gas 
metallicity and density are low ( Gil de Paz et al.|2007 ). On 


the other hand, judging from the observed UV emission, star 
formation might be expected to be supported for at least a 
few billion years. Yet the efhciency of transforming gas into 


a young stellar population is expected to be low (Bigiel et al. 


2010), with a time scale ~ 50 times longer than in normal 


disks. Note that cluster counts in the outskirts of M 83 ap¬ 
pear to be consistent with model predictions based on a 
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standard IMF and cluster aging effects; therefore even low- 
mass newborn stellar clusters occasionally have 

O stars ( Koda et al.|2012 ). 

The thickness of the disc is likely to increase with radius. 
A proper model of this feature and of the dynamical mecha¬ 
nisms involved in star formation and spiral structure noted 
above would also depend on the characteristics (amount and 
spatial distribution) of the dark matter present, which are 
not easy to measure and are obviously of great astrophysical 
interest. 

Thermal instabilities are probably less important, be¬ 
cause of the low metallicity of the gas. Thus, under these 
conditions, to produce star-forming regions a mechanism 
is needed able to compress the gas on scales of about one 
kiloparsec. This suggests that a hierarchy of perturbations 
should be operating, allowing gas to form stars at small 
scales. The perturbation on the largest scale is the one as¬ 
sociated with grand-design spiral patterns. 

On the one hand, natural mechanisms for the genera¬ 
tion of structures beyond the optical radius are those re¬ 


lated to gas accretion (Dekel et al. 2009), streams from 
gas-rich companions in close passages ( Dobbs et al.|[20T0 ), 
and, in general, tidal interactions with satellites or dwarf 
galaxies ( Bullock &: Johnston||2005 ). In these scenarios, the 
role of dark matter is not clear. On the other hand, it is 
generally believed that, in the presence of a triaxial (i.e.. 


non-axisymmetric) distribution of dark matter (Khoperskov 
et al. 2012| Valenzuela et al. 2014), disc-halo interactions 


could generate spiral structure. However, for the majority 
of galaxies our knowledge of the three-dimensional distribu¬ 
tion of dark matter is quite poor. 

Hydrogen is mostly neutral in the outermost regions, so 
that the relatively low photoionization level expected makes 
it likely that large-scale magnetic fields should not play an 
important role in structure and star formation. Although 
it is generally believed that gravitational instabilities are 
the main driver for spiral structure formation, the collective 
processes occurring in the outer gaseous medium of galaxies 
remain largely unexplored. 

The study of the global large-scale spiral structure has 
a long history. The commonly accepted picture is that spi¬ 
ral structure is the manifestation of density waves in the 
galactic disc ( Lin Shu||1964 ), for which gas and stars co¬ 
operate collectively ( Lin V Shu|1966 ). Currently, the picture 
that the grand-design spiral patterns are associated with few 
self-excited global spiral modes is the scenario that has been 
worked out in greatest quantitative detail (see [Bertin Lin 


1996| and references therein), supported by a number of suc¬ 
cessful observational tests. Convincing tests from realistic 
numerical simulations remain difficult to obtain, because of 
the complexity of the physical phenomena involved, ranging 
from the role of resonances in the collisionless stellar com¬ 
ponent to the destabilizing and self-regulating role of the 
dissipative cold interstellar medium (Elmegreen & Thomas- 


son|1993 Baba et M^|2009 Fujii et al.||2Qll ). 

In this paper we extend the study of spiral patterns 
outside the optically-bright stellar disc presented by [Bertin] 
& Amorisco (2010). We perform 3D hydrodynamical non¬ 


linear simulations in simple galaxy models outside a central 
disc; we assume that in the central disc spiral structure is 
dominated by one or few modes, which act as a central “en¬ 
gine” for what is observed in the outer parts. We then study 


the properties of spiral structure in the outer gaseous disc 
by varying a set of parameters that characterize the relevant 
perturbations. If the conditions in the galactic disc favour 
the leakage of small-amplitude (^ 10%) non-axisymmetric 
density waves through the outer Lindblad resonance (OLR), 
we find that indeed these perturbations can give rise to 
large-scale prominent spirals covering a wide radial range 
beyond the optical radius of the galaxy. Note that outside 
OLR the gas can support density-wave propagation even 
when its effective velocity dispersion is above the condition 
of marginal stability with respect to axisymmetric perturba¬ 
tions (or, correspondingly, the column density is below the 
related critical density). 

The paper is organized as follows. Section 2 describes 
the adopted basic model and numerical approach. In Sect. 3, 
we present the main results for various models. In particu¬ 
lar, we address models in which the inner disc is dominated 
by a single mode and models with two or three important 
modes, incorporated as a boundary condition at an inner an¬ 
nular region in the outer parts of the stellar disc. The impact 
of small-scales inhomogeneities and subgrid physics are also 
discussed briefly. Our nonlinear simulations are compared 
with the results of the linear theory. Section 4 provides dis¬ 
cussion and conclusions. 


2 MODEL 


The pure hydrodynamical approach based on the TVD 
MUSCL (Total Variation Diminishing Multi Upstream 
Scheme for Conservation Laws) scheme is used for numeri¬ 
cal simulations. The computational technique was described 
by Khoperskov et al. ( 2014| . Here we present a set of 3D sim¬ 
ulations on a uniform 2048 x 2048 x 128 Cartesian grid. For 
a fiducial model, the computational box is 144 x 144 x 9 kpc 
with a spatial resolution of about 70 pc. A study of reso¬ 
lution effects has been performed with cell size of 200 and 
35 pc (see Sect. [3?7| . The general set up of the simulations 
is illustrated in Fig. 


2.1 Basic state 


Below we deal with the evolution of a gaseous disc embedded 
in the fixed external potential of an approximately isother¬ 


mal dark matter halo (model introduced by Burkert 
1995 ) combined with the potential of a stellar disc 4>d of 
the Miyamoto-Nagai form ( [Miyamoto fc Nagai|1975 ). In di¬ 
mensional units we assume that h = 3 kpc is the stellar disc 
exponential scale length. 

The parameters of the external potential were cho¬ 
sen to support the flat rotation curve of the gaseous 
disc (see Fig.[^. For the majority of nearby galaxies, the gas 
distribution follows an exponential profile within the optical 
radius ropt ( [Bigiel Blitz|20T2 ). At far away distances, gas 
density profiles are not so well known. Nevertheless, it is nat¬ 
ural to assume a Eg oc 1/r profile of gas beyond (1 — 2)ropt. 
For simplicity, we adopt as initial gas density distribution 
the Eg oc 1/r profile also for the inner part of the disc, 
which is kept fixed in our simulations. The initial stellar 
disc surface density and gaseous surface density profiles are 
shown in Fig. [^ 

For the value of the gas velocity dispersion we adopt 
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Figure 1. General scheme of the computational model. A fixed 
area of initial conditions inside the optical radius, defined here as 
^opt = 6h, is shown as a red disc. The thin, light circular annulus 
(cyan line) is the site where the perturbation is imposed, which is 
the inner boundary for the evolving domain. As shown on the top 
right corner, the cyan line corresponds to the two computational 
grid cells where the initial conditions are designed to mimic the 
presence of an outgoing density wave signal coming from the inner 
regions where coherent spiral structure is present. The blue area 
outside the red disc is the actual computational domain from 6h 
out to 24h, where the hydrodynamical quantities evolve. 



Figure 2. Adopted rotation curve V{r). 


that obtained from the marginal-stability condition for a 
non-self-gravitating layer of finite thickness, as is expected in 
the galactic outskirts (see Eqs. (15.21), (15.22), and Fig. 15.4 
in |Bertin|2014l l: 


Q = Qmax = 0.425 . 


( 1 ) 


Thus, using the epicyclic frequency k, and from our ini- 



Figure 3. Initial surface density distributions of gas (solid line) 
and stars (dashed line). 


tial conditions we have the relation for the radial velocity 
dispersion: 


c = 0A257rG'Eg/K , 


( 2 ) 


which, outside the circle r = 6/i, is approximately constant 
c ~ 3.4 km s“^. Such value is consistent with the observa¬ 
tional data, which suggest a velocity dispersion of the gas 
clouds in the range 1 — 10 km s“ 

Elmegreen Scalo||2004 ). 


(Stark & Brand 


1989 


The gaseous disc thickness is set by the condition of 
vertical hydrostatic equilibrium: 


1^ 

p dz dz ’ 


( 3 ) 


where the gas volume density p and pressure p are connected 
by the equation of state p = pc^ and the total gravitational 
potential $ = takes into account the potential 

of the gas which is the solution of the Poisson equation: 


l_a 

r dr 


d^c 


dr 


d^^a 


+^=47rGp. 




( 4 ) 


Eqs (|^ - (|^ determine the equilibrium vertical distribu¬ 
tion of the gas. We define its vertical scale height zo{r) by 
minimizing the difference F(zo): 


F{zo) = Ipocosh ^{zlzo) - p{r,z) \ . 


( 5 ) 


As the vertical gravity decreases with radius then the 
equilibrium disc thickness increases. Figure shows such 
disc flaring from the solution of Eqs. § - @ for our 
initial parameters. We also compare our solution with 
the gas thickness profile supported by the gas self-gravity 
alone (see Eq. (14.12) in Bertin|2014 ): 

ttGE, 


— = 0.18 
r rhi- 


where the coefficient 0.18 is obtained from the marginal sta¬ 
bility condition (see Fig. 2 in Bert in Amorisco|2010 ). 


2.2 The imposed perturbations 

We now introduce the properties of the perturbation im¬ 
posed at the inner boundary (see Fig.[^. Basically, we con¬ 
sider the perturbations of the hydrodynamical quantities X 
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Figure 4. The solid line is the gaseous disc scale height obtained 
from the vertical equilibrium according to Eq. The dashed line 
is the disc thickness according to the marginal-stability condition 
of a fully self-gravitating slab only (Eq. [^. 


at the inner boundary (cyan thin circular annulus in Fig.[^ 
to be proportional to cos {mO — cut) , where m is the mode 
azimuthal number, Qp = Lu/m is the angular speed of the 
spiral density perturbation, t is the time, 0 is the angular co¬ 
ordinate. We assume that the relative amplitude of the den¬ 
sity wave is equal to Aq (mean values are shown in Table 
and then relative amplitudes of the perturbation for all other 
quantities can be found straightforwardly from Eqs. (7), (8), 
(9), and (13) for the short-trailing wave-branch in [Bertin] 
Amorisco| ( |2010 ). In the calculation, the perturbation at 
the inner boundary is applied only for the hydrodynamical 
evolution of the outer disc. Thus the gas distribution is kept 
“frozen” and axisymmetric inside the disc dehned by the cir¬ 
cle Topt = 6/i. For any given mode considered in our study, 
the value of Qp sets an outer Lindblad resonance inside the 
inner boundary at ropt = 6/i. Thus the perturbations that 
we consider propagate outward (see Fig. [^. 

It is believed that several spiral modes generally coexist 
given disc model ( Bertin et al.||197'7| [Korchagin et al. 


2000). Then we consider not only a model where a single 


dominant mode is present, but also models that include a 
superposition of more than one mode. In this case the den¬ 
sity perturbations can be written in the following form: 


Fi = ^ Ao^i cos {rriiO — ujA + 5i) 


( 7 ) 


where Ap^i are the relative amplitudes of the density per¬ 
turbation at the inner boundary, 5i are the initial phases of 
perturbations and uji = Qp^i/rrii. For the simplicity, in the 
following we take (5^ = 0. Much like for the case of a single 
mode perturbation model, amplitudes for pressure and ve¬ 
locities are calculated according to the expressions provided 


by Bertin & Amorisco (2010). We consider three models 


with multi-mode boundary conditions (see detailed param¬ 
eters in Table [^. In particular, we take a case that includes 
the superposition of a one-armed and a two-armed pertur¬ 
bation (model HI), a case with a pair of two-armed pat¬ 
terns (model FI), and the case made of three modes, each 
with a different number of arms (model JI). The amplitudes 
and the corresponding pattern speeds of the perturbations 


Figure 5. Kinematics of the gaseous disc. The thick solid line 
is the angular velocity Q = V{r)/r, the thin solid lines represent 
H±Av/ 3, the dashed lines Q^k,/2 and the dash-dotted line is 
The horizontal thin solid lines correspond to the pattern speeds 
of the modes considered in our simulations: from top to bottom, 
50 km s“^ kpc“^(for the m = 3 mode), 40 km s“^ kpc“^(for 
one m = 2 mode), 35 km s“^ kpc“^(for a second m = 2 mode), 
30 km s“^ kpc“^(for the m = 1 mode). The mode parameters are 
summarized in Tabled All the selected modes have OLR inside 
the inner boundary of the computational domain (see also Fig.[^. 


for these models are chosen so as to be in qualitative agree¬ 
ment with the linear theory of global spiral modes. Of course, 
other combinations of perturbation parameters could be rea¬ 
sonable: our models should only be considered as a simple 
representation of typical cases broadly consistent with the 
modal density wave theory (seejBertin et al.|I977||Bertin & 
Lin|1996| |. 


To avoid an initial exaggerated kick on the disc, we let 
the perturbation amplitudes grow from vanishingly small to 
hnite values according to a linear law during one typical 
rotation period Ti until they reach the chosen value Aq. 
From Fig. our time scale is the rotation period Ti 
0.5 Gyr at 6h = 18 kpc. 

Initially we set up the dynamical equilibrium of the 
gaseous disc within the entire computational domain, that 
is, in both the red and the blue areas of Fig. Then we 
ignore the complex self-consistent evolution of the stellar- 
gaseous disc within the optical radius. Namely, we keep the 
central part of the computational domain r < ropt = 6/i as 
“frozen” (see Fig. [^. This dehnes an inner boundary layer 
at ropt = 6/i for the live gaseous disc outside. 

It is believed that the rotation of the galaxy outside 
the bright optical disc is mainly supported by the gravi¬ 
tational potential of the dark matter, which in turn is ex¬ 
pected to have signihcant substructures ( Moore et al.|I999 ). 
In fact, numerical simulations of galaxy and structure forma¬ 
tion in the cosmological context predict that galactic dark 
matter haloes should contain a population of so-called sub¬ 


haloes (Giocoli et al. 


2008 


Gao et al. 2004). The relative 


motions of these substructures should induce local time- 
dependent variations of the gravitational held. This feature 
might have a strong impact at the periphery, where the bary- 
onic matter density is small. To consider these effects we 
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Table 1. Parameters of different runs, where m is the perturbation number of arms, Aq is the relative amplitude of the adopted density 
perturbation, Qp is the pattern speed, Tco is the corotation radius, roLR is the outer Lindblad resonance position, h = 3 kpc is the 
exponential scale length of the stellar disc. When several parameters are present in a given column the inner boundary is perturbed with 
more than one mode. B1 is the reference model. 


Run 

m 

Ao 

Up 

km kpc“^ 

r CO 

h 

^OLR 

h 

Additional 

feature 

EI 

I 

0.05 

30 

2.4 

5.6 

- 

BI 

2 

O.I 

40 

1.8 

3 

- 

B2 

2 

O.I 

40 

1.8 

3 

clumpy gas distribution 

B4 

2 

O.I 

40 

1.8 

3 

potential perturbation 

B5 

2 

O.I 

40 

1.8 

3 

subgrid cooling 

B7 

2 

O.I 

40 

1.8 

3 

higher velocity dispersion c = 5 km s“^ 

KI 

3 

0.15 

50 

1.4 

2.1 

- 

HI 

1/2 

0.05 / O.I 

30 / 40 

2.4/1.8 

5.6 / 3 

- 

FI 

2/2 

0.07 / O.I 

35 / 40 

2/1.8 

3.4 / 3 

- 

JI 

\ 12 13 

0.05 / O.I / 0.15 

30 / 40 / 50 

2.4/I.8/I.4 

5.6 / 3 / 2.1 

- 


performed also simulations with a time-dependent inhomo¬ 
geneous gravitational field. 

The simulation of the dark matter dynamics at the out¬ 
skirts would require high spatial and mass resolution. To 
avoid the related technical problems, in our numerical model 
we add to the initial conditions adopted for B1 a random 
perturbation of the halo gravitational potential. This de¬ 
fines model B4. That is, in the B4 model we recalculate the 
potential according to the rule: 


$/i(r, t) = ^h{r, z, 0) [1 + a{t)] 


( 8 ) 


where a{t) is a random number in the range [—0.1; 0.1] which 
varies at each time step and is unique for cells with size 
200^ pc^. Thus the total mass of the halo is conserved but 
a 10% time-dependent perturbation of the gravitational po¬ 
tential is introduced. 

The galactic gaseous component is a cloudy medium. 
Within the optical disc a significant fraction of the gas mass 
is concentrated in giant molecular clouds. Outside the opti¬ 
cal disc the cold gaseous phase is likely to be concentrated 
in neutral hydrogen clouds (but see the picture explored 
by Pfenniger et al. ( 1994| )); apparently, these clouds are not 
forming stars in large amounts. These arguments suggest 
that we should consider an inhomogeneous gas distribution 
in our simulations. We designed model B2 so as to include 
a random perturbation of the gas density distribution with 
relative amplitude of 10%. The velocity field of the clouds 
was perturbed according to the mean velocity dispersion of 
about 3.4 km s“^. As was mentioned earlier, this value ap¬ 
pears to be realistic. 

It is generally thought that the gas cloudy medium is 
collisional. To take into account this fact, we calculate an ef¬ 
fective “cooling rate” of the cloudy medium associated with 
inelastic collisions. We assume that in each computational 
cell there is a subgrid population of clouds that can lose 
energy as a result of collisions. Obviously the cooling rate 
depends on the gas density n and cloud velocity dispersion 
at the cell. In calculations we applied the cloud collision rate 


based on the model by Ricotti & Ferrara (2002). In the nu¬ 


merical scheme the cooling rate was used as a source term 
in the energy conservation equation. It is implemented with 
the standard technique used for the radiative cooling ap¬ 


proximations widely adopted in simulations of galaxies and 
ISM. 

In the next sections we consider the results of the dy¬ 
namical simulations. First we discuss the reference case of 
a model with a single m = 2 mode. We pay attention to 
the spiral morphology, its time-dependent evolution, and we 
check to what extent its behaviour agrees with the linear the¬ 
ory (Sects. [3J^ 3.3). Next we consider a more complex and 
realistic situation by studying the case of a clumpy gas distri¬ 
bution and of an inhomogeneous gravitational potential. A 
model with subgrid energy dissipation resulting from inelas¬ 
tic HI cloud collisions is also described in Sect. |3.5[ Finally 
we study cases with multi-mode perturbations (Sect. [3^. 


3 RESULTS 

3.1 Single-mode perturbation 

We now describe the results of hydrodynamical simulations 
of the gaseous disc evolution with a single-mode m — 2 
perturbation imposed at the inner (^opt) boundary. The 
values of the adopted parameters are shown in Table 
In the following, the total density of the gas is denoted 
by — (E^) + El, where (E^) is the azimuth-averaged 
(radius-dependent) profile of E^ and Ei represents the spi¬ 
ral density wave. In general we will refer to E^ in units of 
Mqpc“^; on other occasions, such as in Fig. 6, we will refer 
to the (dimensionless) relative density perturbation defined 
as El/(Eg,). 

First we discuss the reference model BI with azimuthal 
number m = 2, pattern speed Qp = 40 km s“^ kpc“^, and 
relative density amplitude Aq — O.I. In Fig. the evolu¬ 
tion of the surface density perturbation is shown. The two- 
arm trailing spiral structure imposed at the inner boundary 
moves outwards and its amplitude increases in the outer 
parts. At time ~ Ti after the beginning of the simulation, a 
quasi-stationary structure sets in. 

In the inner regions, inside ~ 2ropt, the spiral density 
waves are regular and characterized by low relative ampli¬ 
tude ^ O.I — 0.3. In the outer parts, the wave amplitude 
increases rapidly; the shape of the density perturbation de¬ 
parts from being sinusoidal and the density perturbation be- 
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time [t/T^] 

Figure 6. Evolution of the relative surface density perturbation (Ei/(E^}) in model B1 at times t = 1; 1.1; 1.2; 2; 2.7; 3 in units of 
500 Myr. The plots are drawn in the inertial, nonrotating frame of reference. Black circles are drawn at radii 6h, 12h, 18h, and 24h. 
The bottom right panel shows the evolution of the amplitude of the m = 2 Fourier component of the relative surface density perturbation 
at different radii, marked by lines with thickness decreasing with increasing radius. 


10 rr-r-r- 


10 r 


—1— \ —I—I—I— \ —I—I—I—I— I —I— I —I—I—I—I—I—I—I—I—I—r— 


instability is clearly seen in the evolution of the amplitude 
of the m = 2 Fourier component in model B1 (see Fig. |^. 
In general, this amplitude remains approximately constant 



r=7h ' 

at all radii after the initial transient, that is, at times longer 
than t — 1.8 — 2Ti. 

The origin of the shock instability is likely to be similar 
to that of the the wiggle instability investigated numeri- 


.1 ^ 


i ' 

r=l^ 

cally by several authors (e.g., Wada & Koda 2004| Dobbs 

1_1_1_1_1_1_1_1_1_L 


j r=16h_ 

r=20h' 

& Bonnell||2006| Kim et al.||20I4). In other words, Kelvin- 
Helmholtz instability is likely to be one important mecha¬ 
nism for the formation of spurs in the vicinity of the spiral 
shocks. The numerical simulations by Wada h Koda (2004) 

r=23h' 

_i_1_1_1_1_1_1_1_1_1_1_1_1_1_1_ 

suggest that tightly-wound spiral shocks should be relatively 
stable, compared with the case of open spirals. However, in 


100 200 
e fdegl 


300 


Figure 7. Total surface density profiles for the gas (E^), in units 
of Mqpc“^, are shown along the azimuthal coordinate at given 
radii for model B1 at t = 2.5Ti. The line thickness decreases with 
increasing radius. 


comes asymmetric (see Fig.0. Beyond ^2 — 2.5ropt rather 
narrow shocks form. In fact, the Mach number \y (r)—rQp]/c 
increases linearly with radius because c ^ 3A km s“^and 
V{r) 210 km s“^are approximately constant (see Eq. ^). 
Thus the spiral structure is characterized by supersonic mo¬ 
tion across the disc. 

In the outermost regions, where the relative amplitude 
of the spiral density perturbation Ei/ (E^) reaches the values 
of ~ 1 — 2, then the shocks become unstable and small-scale 
spurs appear. It is likely that the instability of spiral shocks 
makes the perturbation saturate at finite amplitudes. This 


our simulations the large amplitude of the shock associated 
with the observed instability occurs in the outermost parts 
of the disc (r > 12h ~ 2ropt), where the pitch angle is 
smaller than 5° (the radial profile of the pitch angle i{r) is 
illustrated in Fig.[^. In turn, in the inner parts of the com¬ 
putational domain (r < 12h) the pitch angle is close to 10°; 
but there the amplitude of the gas perturbation is relatively 
small so that the density waves are in the linear regime (see 
also the caption to Fig. 0 and the strong shock instability 
is suppressed. It should be noted that in our simulations the 
spatial size of spurs is about one kiloparsec, which is com¬ 
parable to the size of some small-scale structures inside the 
optical radius of nearby galaxies. 


3.2 Direct measurements of pitch angles and 
pattern speeds 

To derive the radial profile of the pitch angle of the spiral 
arms observed in our simulations, we have used the method 
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Figure 8. Directly measured pattern speed at different radii for 
model Bl. 


described by Davis et al. (2012). We consider the Fourier 


analysis of the surface density distribution of the perturba¬ 
tion El as a function of azimuth 0 and the Fourier transform 
in the logarithmic radial coordinate u = logr (this is often 
referred to as a decomposition in logarithmic spirals). The 
integrations are performed over narrow annular rings in the 
computational domain, that is: 


1 nri + l n 

A{p,m,ri) = TT / / 

^0 J Ur- j — TT 


0)e 


— i{m9-\-pu) 


dOdu , (9) 


r'^24h 

where Go = / Ei(r^, 0)d0du is a suitable normaliza- 

tion constant. Thus A{p,m,ri) represents the contribution 
of the m-armed logarithmic spiral component, with pitch 
angle i = arctan (—m/p) at given radius n. By considering 
the value pmax at which the quantity A{p,m,ri) attains its 
maximum at given m and ri we can thus reconstruct the 
pitch angle radial prohle for our simulated spiral structures. 
The error on the pitch angle value ii can be found from the 
spatial variation of the quantity in the range ri+ 1 / 2 ]. 

In the narrow annuli at radius ri considered in the 
method just described, we can also measure the pattern 
speed associated with the spiral structure present. To make 
such measurement, we proceed as follows. A phase angle for 
given m and p dehning the orientation of the spiral pattern 
at radius can be calculated as 


T = arctan 


Im(A) 
Re(A) ’ 


( 10 ) 


where Im(A) and Re(A) are the imaginary and the real part 
of A{pmax,rn), respectively. Then a local value of the speed 
of the pattern with given m can be determined as 

1 

= ( 11 ) 
m at 

The error on the value of the pattern speed depends on both 
the spatial variation of the derivative in Eq. (11) and on 


small time-dependent variation of the quantity that is cal¬ 
culated. When a single mode is imposed at the inner bound¬ 
ary, the procedure indeed gives back the value of the pattern 
speed of the imposed perturbation. In Fig. calculated at 



r/h 


Figure 9. The relative amplitude of the density wave in model 
Bl is shown by squares. Circles represent the radial velocity per¬ 
turbation relative to the rotation velocity of the basic state for 
the same model. The radial profile of the pitch angle in model Bl 
is shown by error bars; error bars are associated with the averag¬ 
ing of the pitch angle in rings of finite radial size (see Eq. [^. T he 
solution from the linear theory by [Bertin &: Amorisco| l |2010^ is 
shown by lines. 


t ^ 2Ti for the Bl model, the measured pattern speed is 
shown to be constant with radius and consistent with that 
of the single-mode perturbation imposed at the inner bound¬ 
ary. 


3.3 Comparison with the linear theory 


From the linear theory of density waves [Bertin fc Amorisc^ 
(2010) obtained the expressions for radial velocity, surface 
density, and pitch angle of short-trailing density waves in the 


galactic outer regions (e.g., see Eqs. (9) and (13) in Bertin 
Amorisco|2010 ). In this section we compare the results of 


our dynamical simulations with the predictions of the linear 
analysis. 

Eor the frame at t = 3Ti from the reference Bl simula¬ 
tion (see bottom row in Eig.|^ we calculate the radial prohle 
of surface density perturbation, radial velocity perturbation, 
and pitch angle (see Eig. |^. The perturbation amplitudes 
of surface density and radial velocity increase with radius, 
whereas the pitch angle of the pattern decreases. This quali¬ 
tative behaviour agrees well with the expectations of the lin¬ 
ear theory (see Eig. 5 in Bertin Amorisco|2QlQ ). Quantita¬ 
tively, the good agreement between linear theory and simula¬ 
tions applies to a rather wide radial range 6h < r < 10 —12/i. 
At larger radii (r > 2ropt), the simulations exhibit a strongly 
nonlinear behaviour, because the relative perturbation am¬ 
plitudes attain very high values, up to 2 — 3. 


3.4 Model with higher velocity dispersion 

In this section we consider a different model (B7) charac¬ 
terized by a velocity dispersion higher than that adopted in 
the reference Bl model, well above the value required by 
the condition of marginal axisymmetric stability. We thus 
take c = 5 km s“^. Interestingly, the dynamical evolution of 
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this hotter system basically follows the same picture as de¬ 
scribed in Sect. |3.1| Of course some morphological changes 
are expected and indeed found in the simulations. 

In Fig. we show the relative surface density distri¬ 
bution established at time 2.5 Ti. Because of the higher gas 
velocity dis- persion, in the B7 model we observe the ex¬ 
citation of a significantly more open spiral structure (with 
respect to the reference B1 model). In the regions close to 
the inner boundary, the pitch angle is 20° (to be com¬ 
pared to the value of ~ 10° found in Bl). Because a larger 
pitch angle of the pattern provides better conditions for the 
shear instability of the spiral shocks, a more perturbed mor¬ 
phology of the pattern and the presence of prominent spurs 
and feathers along the spiral arms are expected and found 
in the simulations. 


3.5 More realistic models 


Our reference model Bl demonstrates the possibility of reg¬ 
ular and sharp spiral patterns, of the type that is observed in 
some deep HI images (e.g., in NGC 1512), but physically it is 
exceedingly simple. More realistic models should be devised. 
In particular, we have checked how different small-scale pro¬ 
cesses, which can be implemented at the subgrid level in our 
simulations, affect the morphology of the observed spiral 
structure. 

Here we compare four types of models with the same 
imposed perturbation at the inner boundary: (i) the refer¬ 
ence Bl model, which was described previously; (ii) the B2 
model, which is based on a clumpy gas distribution and a 
smooth potential; (iii) the B4 model, which is based on a 
smooth gas distribution and a halo potential perturbed by 
clumps of dark matter; and (iv) the B5 model, which is 
the same as Bl, but takes subgrid cooling into account. Re¬ 
sults of simulations for these models are shown in Fig. 
where the gas surface density perturbation Ei is illustrated 
at t ~ 2.5Ti. The basic conclusion is that the results shown 
for the Bl model are robust. 

In fact, the large-scale morphology (grand design) is 
similar for all models. However, various additional features 
are found to characterize the small-scale morphology. For 
the B2 and B4 models the new small-scale features basically 
cover the entire disc. In contrast, in the B5 model, with 
subgrid cooling, the effects are most evident in the denser 
regions at the edge of the spiral shocks, making the arms 
less smooth; the intensity and spatial scales of the spurs 
and feathers that are observed in the simulations suggest 
that the B5 model involves the wiggle instability known to 
affect a multi-phase inhomogeneous interstellar medium in 


the presence of spiral shocks on the galactic scale (Wada 

200 ^. 


several non-axisymmetric modes is imposed. The HI model 
is based on the combination of an m = 1 mode rotat¬ 
ing at 30 km s“^ kpc“^and an m = 2 mode rotating at 
40 km s“^ kpc“^. The FI model considers the superposition 
of a pair of m = 2 modes with different amplitudes (^o,i = 
0.05, Ao ,2 = 0.1) and pattern speeds (30 km s“^ kpc“^, 
40 km s“^ kpc“^). The J1 model studies the case in which 
three modes with different m numbers are present (m = 1, 
m = 2, and m = 3). A description of the adopted param¬ 
eters is given in Table Figure illustrates the surface 
density perturbation maps for the single-mode cases, that 
is, the standard Bl model (m = 2), the El model {m = 1), 
and the K1 model (m = 3), and for the models with several 
modes (HI, FI, Jl). 

The evolution of the disc in the case in which several 
modes are imposed at the inner boundary is rather simi¬ 
lar to that of the single-mode case. The spiral morphology 
of discs where a dominant m = 3 perturbation is applied 
tend to exhibit a prominent three-armed structure during 
the simulation. The presence of an m = 1 imposed pertur¬ 
bation is generally associated with some lopsidedness. The 
models are characterized by a time-dependent evolution of 
spiral structure, although the overall observed patterns do 
not appear to vary very significantly in time. Of course, the 
evolution is associated with the superposition of the modes, 
rotating with different angular speeds. The qualitative be¬ 
haviour remains generally similar to that of the Bl model, 
with nonlinear behaviour setting in at large radii and some 
wiggle instability occurring when subgrid cooling is incor¬ 
porated. 


3.7 Resolution study 


All the models described so far have the same spatial reso¬ 
lution, that is, a cell size of 70 pc in physical units. In order 
to check whether spatial resolution affects our general re¬ 
sults, we have performed simulations with lower (200 pc) and 
higher (35 pc) cell linear size. Figure 


13 


shows the results of 
these simulations. The general grand-design spiral structure 
is basically unchanged, as expected. Obviously, the narrow 
shock is smoother in the low resolution simulation; higher 
spatial resolution let us resolve instabilities growing on very 
small scales related to shear flows behind the shock (see also 
Wada fc Koda|2004 |. 


In conclusion, our reference 70 pc resolution is likely to 
be sufficient for studies of the global spiral structure and 
the detection of the shock instability. However, as to issues 
related to local gravitational instabilities and star formation 
processes in the outermost gaseous discs in galaxies, more 
detailed simulations with higher spatial resolution are de¬ 
sired. 


3.6 Simulations with more than one mode 
imposed at the inner boundary 

So far we have described models on which a single mode is 
imposed at the inner boundary. In real grand-design spiral 
galaxies, it is natural to expect that the large-scale mor¬ 
phology is dominated by the superposition of few spiral 
modes, each characterized by its own amplitude and pat¬ 
tern speed. We thus investigate the properties of simula¬ 
tions in which at the inner boundary a superposition of 


4 DISCUSSION AND CONCLUSIONS 

In this paper we have presented a set of 3D hydrodynam- 
ical simulations that describe the establishment of large- 
scale regular spiral structure in the outermost gaseous disc 
of spiral galaxies. The general scenario considers a galaxy 
in which spiral modes are excited inside the bright opti¬ 
cal disc through the transfer of angular momentum to the 
outer regions by means of short trailing density waves. In the 
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Figure 10. Relative surface density perturbation (Ei/(E^)) for different models, with a single m = 2 perturbation imposed at the inner 
boundary, at time t = 2.5Ti. From left to right: top row, B1 (the reference model) and B7 (with higher gas velocity dispersion); bottom 
row, B2 (with a clumpy initial gas distribution), B4 (with a clumpy dark halo), B5 (with subgrid cooling). For a detailed parameter 
description, see Table [T] Black circles mark different radii, at 6/i, 12h, 18h, and 24h. 


gas, such outgoing density waves can leak through the outer 
Lindblad resonance and propagate outwards even when the 
gas effective velocity dispersion is above the condition of 
marginal stability with respect to axisymmetric perturba¬ 
tions (or, correspondingly, the column density is below the 
related critical density). In the simulations of the outermost 
gaseous disc, these outgoing waves are imposed as a station¬ 
ary disturbance at the inner boundary. We have investigated 
the case in which galaxies and the corresponding boundary 
conditions are dominated by a single mode and, separately, 
the case in which more than one important mode is present. 
The results that we have obtained can be summarized as 
follows: 

(i) The simulations are run as studies of a time-evolving 
situation in which the inner boundary acts as a source of 
density waves. After a relatively rapid initial transient, a 
quasi-stationary spiral structure is established over the en¬ 
tire outer disc, well outside the bright optical disc. 

(ii) The simulations exhibit very good quantitative 
predictions of the linear theory by |Bertin| 

out to r ~ 1.5ropt- At larger radii the 
amplitude of spiral structure increases beyond the reach of 
the linear theory and then saturates. In this sense, our model 
is reminiscent of nonlinear tsunami-like waves. Correspond¬ 
ingly, spiral shocks form in the outermost regions, as a result 


agreement with the 
& Amorisco (2010 


of the supersonic motion of the patterns through the gaseous 
medium. Spiral shocks tend to be Kelvin-Helmholtz unsta¬ 
ble in the post-shock regions. 


(iii) The simulations suggest that the outer spiral struc¬ 
ture may be associated with significant star formation. As 
in a galactic tsunami, small amplitude perturbations be¬ 
come stronger and stronger at large radii, so that the shocks 
formed might trigger star formation events; some indications 
of star formation are indeed noted in UV observations. We 
note that the process studied in this paper might explain 
an outer UV star-forming ring in isolated galaxies even in 


the absence of an impact by an external object (Ilyina et al. 


2014). In our simulations, we did not investigate the star 


formation processes and the issue of the expected UV flux 
in great detail. The main reason for this is that to carry 
out a proper investigation of these processes would require 
a deep study of the conditions of the gaseous medium in 
the galactic outskirts, which are at present not well con¬ 
strained by the observations; furthermore, we should have 
dealt with issues related to the resulting IMF, which are 
even less known. In this respect, we think that producing 
from the simulations synthetic UV spectra to be compared 
with the observations would be premature. In turn, in our 
simulations we focused on the larger-scale dynamical aspects 
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Figure 11. Relative surface density perturbation (Ei/{E^}) for models with different conditions imposed at the inner boundary, at 
time t = 2.5Ti. Top row: El (m = 1), B1 (m = 2), K1 (m = 3). Bottom row: HI (superposition of an m = 1 and an m = 2 mode), 
FI (superposition of a pair of m = 2 modes), J1 (superposition of an m = 1, an m = 2, and an m = 3 mode) (see Table[^. Black circles 
mark different radii, at 6h, 12h, 18h, and 24h. 




Figure 12. Evolution of the relative surface density perturbation (Ei/(E^}) in model J1 with three modes imposed at the inner boundary, 
at times t = 1.1; 1.4; 1.7; 2; 2.3; 2.6; 2.9; 3.2 in units of 500 Myr. Black circles are drawn at radii 6h, 12h, 18h, and 24h. The plots are 
drawn in the inertial, nonrotationg frame of reference. 
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Figure 13. Relative surface density perturbation (Si/ (S^)) at t = 2.5Ti for the B1 model simulated with varied spatial resolution, that 
is, different linear cell size: left — 200 pc (low resolution), center — 70 pc (reference case), right — 35 pc (high resolution). Black circles 
mark different radii, at 6/i, 12h, 18/i, and 24/i. 




of the processes involved during the propagation of density 
waves in the outermost gas disc. 


(iv) The simulations of more realistic models, includ¬ 
ing a variety of physical factors, exhibit a generally simi¬ 
lar behaviour in relation to the large-scale spiral structure. 
This suggests that the results obtained are rather robust and 
that indeed prominent large-scale spiral patterns should be 
a natural feature of galaxies with a gaseous disc extending 
beyond the optical radius. The simulations by |Bush et al.| 
( 200^ (see also [Bush et al.||201Q| are interesting and ex¬ 
hibit some features in common with the results of our paper 
(in particular, the frequent finding of organized compres¬ 
sion regions and filamentary structures). However, we wish 
to note that the study by Bush et al.| ( [200^ : (1) focuses on a 
’’fiducial star formation law”; (2) pays little attention to the 
relation between spiral structure in the bright optical disc 
and spiral structure in the outer gaseous disc; (3) appears to 
support the picture that star formation is expected only if 
the gas layer is close to conditions of local Jeans instability; 
(4) appears to depend on the addition of an extended outer 
disc with a constant density; (5) does not discuss the role of 
the thickness of the gaseous layer (which is a crucial factor 
in determining its stability). 


(v) In the picture proposed in this paper, the level of 
structuring of star formation regions in the outermost disc 
should reflect the level of regularity of the spiral structure in 
the bright optical disc. In other words, the spiral structure 
observed in the gas outside the bright optical disc should 
be characterized by well-organized, structured spiral pat¬ 
terns and filamentary structures if the bright optical disc 
is dominated by a grand-design structure. In contrast, the 
outermost spiral arms are expected to be less structured if 
the bright inner disc is dominated by many modes. This is 
a prediction that could be tested by studying a sufficiently 
large sample of spiral galaxies with extended gaseous discs. 
Of course, such a study would go well beyond the scope of 
this paper. In addition, we note that our simulations sug¬ 
gest clearly the possibility of ring-like star-format ion regions, 
because we showed that the pitch angle of spiral patterns 
tend to decrease with radius; evidence of the relation be¬ 
tween these morphological aspects and the underlying flows 


is likely to show up as brighter HI and, possibly, UV emis¬ 
sion. 

Of course, there are a few other interesting issues that 
are beyond the scope of this paper and await further in¬ 
vestigations. One of these issues is the possible use of the 
observed spiral structure in the context of this paper to di¬ 
agnose the amount and distribution of dark matter in the 
outer regions. To this purpose, this study should be further 
validated by detailed comparisons with observations in in¬ 
dividual objects. 
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